Two site self consistent method for front propagation in reaction-diffusion system. 



o 
o 



Niraj Kumar and Goutam Tripathy 
Institute of Physics, Sachivalaya Mary, Bhubaneswar 751005, India 

We study front propagation in the reaction diffusion process A <-> 2A on one dimensional lattice 
with hard core interaction between the particles. We propose a two site self consistent method 
(TSSCM) to make analytic estimates for the front velocity and are in excellent agreement with the 
simulation results for all parameter regimes. We expect that the simplicity of the method will allow 
one to use this technique for estimating the front velocity in other reaction diffusion processes as well. 
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Front propagation in the reaction diffusion system is 
an important field of study in nonequilibrium physics. In 
many natural phenomena we often encounter propagat- 
ing fronts separating different phasesQ . Here, in this pa- 
per, we study the front dynamics in the reaction diffusion 
system A «-> 2A whose mean field description is given by 
the following Fisher-Kolmogorove equation 7] for the lo- 
cal density of A particles p(x, t): ^ — D^^ + kip— k 2 p 2 - 
Here D is the bare diffusion coefficient of A particle 
while ki and k 2 are the rates of creation and annihi- 
lation respectively. This equation arises in the macro- 
scopic description of many processes in natural science 
and serves as a generic model to describe front propaga- 
tion in a system undergoing transition from unstable to 
stable state. The homogeneous steady states of this equa- 
tion are p = f 1 (stable) and p = (unstable). Hence if we 
start with an initial condition where both the states co- 
exist, the the stable phase invades the unstable one with 
speed V as a travelling wave of the form : p = p(x — Vt). 
Here, such fronts are referred to as pulled in which the 
leading edge where p « 1 plays important role in de- 
scribing the front dynamics. Such leading edge analysis 
gives V > Vo — 2\Jk\D. However, for steep enough ini- 
tial condition the minimum velocity Vo is selected and it 
is known that convergence to Vq takes place with a very 

slow power law P3E3 : ^(*) ~ V o + f + 0(t^), where, 
ci < 0. 

In this paper, we study the microscopic realization of 
the process A <-» 2A on a lattice. The important feature 
to be noted here is that the Fisher equation fails to deal 
with internal fluctuations arising due to discrete nature 
of the reacting species, especially when occupancy per 
site N is small HQ. When TV — > oo, the mean field 
results are recovered. 

Motivated by the velocity selection principle, in 
which the leading edge plays an important role, we 
propose a procedure called two-site self consistent 
method(TSSCM). In this method, we describe the front 
dynamics by considering the evolution of occupancy at 
only two sites, the first one is the front site while the 
other one is the site just behind it. In other words, in 
the frame moving with the front we study the evolution 
of occupancy at a site just behind it. By applying a self 
consistent approach (explained later) in this method we 
obtain the analytic estimates for the front velocity which 



are in excellent agreement with the simulation results for 
all parameter regimes. In our simulation we consider a 
one dimensional lattice composed of sites i — 1,2, ..L, 
L being the size of the lattice. Initially, we start with 
left half of the lattice filled with A particles while keep- 
ing the remaining right half empty. Each site can ei- 
ther be empty or occupied by maximum one particle i.e. 
hard core exclusion is taken into account. We update the 
system random sequentially where L microscopic moves 
correspond to one Monte Carlo step(MCS). During each 
update we randomly select a site and the particle at the 
site undergoes the following microscopic moves: (1) The 
particle can diffuse to the neighbouring empty site with 
rate D, (2) The particle can give birth of one particle at 
an empty neighbouring site with rate e , (3) The particle 
gets annihilated with rate W, if the neighbouring site is 
occupied. Due to these microscopic processes the right- 
most A particle, which is identified as a front, moves and 
we arc interested in finding the velocity with which it 
moves, which is given byQ: 



V = e-pt(W-D) 



(1) 



Here, p\ is the density at the site just behind the front 
particle. The above expression for the front velocity can 
be understood by visualising the front particle as a ran- 
dom walker moving with rate e + D in the forward di- 
rection and with rate W pi + D(l — px) in the backward 
direction. 

From Eq.Q, it is obvious that we need to know pi 
in order to predict the front velocity V and there is no 
systematic method to find this p\ exactly. However, sev- 
eral approximations have been proposed, for example, in 
, p\ was taken as the bulk density p — ^ppp. This 
approximation is exact for W = D and shows reasonable 
agreement when | W — D | is small but for larger | W — D\ 
we observe poor agreement, see Fig.l|T]l. In order to make 
better estimate for p\ , especially when we are away from 
the special point W = D, in we proposed a method 
which systematically improves analytic estimate to the 
desired degree of accuracy. Here we write a master equa- 
tion in the frame moving with the front particle and study 
the evolution of one(l — 1), two(Z = 2), three(Z = 3) ..etc 
sites behind the front. For example, for I = 1, we study 
the states {01, 11} and for I = 2 we have the following 
set of states {001,011,101,111}, with the rightmost T' 
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FIG. 1: Simulation results for the density profile behind 
the front for e — 0.25, W = 0.125 for different values of 
D . Inset: density profile for e = 0.05, W = for D = 
0.05,0.15,0.30,0.45 from top to bottom. We note that the 
difference between pi or p2 with bulk density increases with 
increasing magnitude of \W — D\ and for W = D all the sites 
are at the same bulk density. 



FIG. 2: Simulation results for the spatial density correlation 
between pair of successive sites, C(i) =< n(i)n(i + 1) > — < 
n(i) >< n{i + 1) >, for e = 0.25, W = 0.125 for different 
values of D. Inset: C(i) versus i for e = 0.05, W = for 
D = 0.30,0.15,0.05 from top to bottom. Here we note that 
the strength of the correlation increases as \W — D\ increases. 



representing the front. For larger valus of \W — D\, we 
need to study the states corresponding to larger values of 
I in order to get better analytic estimates but writing the 
transition matrix corresponding to the master equation 
becomes more difficult as its size increases as 2 l x 2 l . In 
order to make analytically tractable approximation for p\ 
and hence the front velocity we have presented reduced 
three particle representation. 

In 0], Kerstein proposed a two particle represen- 
tation for the reaction diffusion process A — ► 2A, 
where we study the evolution of following set of states: 
{11,101,1001,10001,....}, with rightmost '1' represent- 
ing the front. The important point that we note here 
is that this set of states is closed under the transition 
for W = while for W ^ it is not closed. This re- 
stricts its applicability for nonzero values of W. Apart 
from this, it uses the ansatz pj — po(l — Po) 1 , where 
Pj is the probability that there are j empty sites be- 
tween the first and second particle. This ansatz clearly 
neglects the spatial density correlation which is, in gen- 
eral, not true, see Fig. J5J). In fact this correlation in- 
creases with \ W — D\ and vanishes for the special point 
W = D. The ansatz also assumes the sites behind the 
front to be at the same density p , which is, of course, 
not true, Fig. JJJ. This leads to a systematic error in 
the analytic estimates which increases gradually with in- 
creasing value of | W — D\. However, one can get better 
and better results by studying the states having larger 
number of particles, for example in [(|, we presented 
mixed representation in which we study the following set 
of states {011, 111, 0101, 1101, 01001, 11001, ...} with the 
rightmost T' representing the front and as expected we 
get better analytic estimates as compared to two parti- 
cle represenation. But writing the master equation gets 



complicated as we increase the number of particles in the 
states. 

Here, in this paper we present a simple analytic esimate 
for p\ and hence the front velocity which are in excellent 
agreement with the simulation results. Here we propose 
a self consistent two-site representation, where we write 
master equation in the frame moving with the front par- 
ticle(as in for the evolution of two states {0A,^L4}. 
Here the rightmost A represents the front particle. In this 
truncated representation these two states make transition 
between each other due to the microscopic processes of 
the model and form a closed set under such transitions. 
For example, OA — > AA if the leading particle gives birth 
to its left empty site with rate e. Diffusion of the lead- 
ing particle to its left changes the state OA — ► AA, pro- 
vided the second site(F — 2) behind the leading particle 
is occupied and remains unchanged if it is empty. If the 
probability of occupancy at F — 2 is denoted as p2 then 
the transition OA — > AA occurs with rate Dp2- Sim- 
ilarly, when the leading particle in the realization AA 
gets annihilated the state of realization changes to OA 
if F — 2 is empty with rate W(l — p-i) while state A A 
remains unchanged if F — 2 is ocuupied. Considering all 
such transitions one can write the master equation for 
the state probabilities Pa a and Pa a as, 

Poa - {2D-Dp 2 + 2W)PAA-(2Dp 2 + 2e + ep 2 )P 0A , 
Paa = {2Dp 2 + 2e + ep 2 )P 0A - {2D- Dp 2 + 2W)P A A- 

(2) 

From Eq. @, it is obvious that we need to know p 2 
in order to find the steady state probabilities Pqa and 
Paa ■ In 5], p 2 was approximated as /5=bulk density 
= e _^ w , which works reasonably well for smaller values 
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for \W — D\ but shows poor agreement with the larger 
values of \W — D\. We note that p 2 , in general, is a 
function of all the parameters (e, D, W) of the model. 
We express the dependence of p 2 on the parameters as 
follows: 



P2 = clPaa + bP\ A 



(3) 



Here, the dependence of p 2 on the parameters is ex- 
pressed implicitly by Paa whose solution is to be sought 
and hence we call the method self consistent. Now using 
the fact that when W = D, p 2 — Paa — P — e /( e + D), 
we rewrite equation© as : 



P2 = (l 



be 



D 



)Paa + bP 



AA 



(4) 



Using this value of p 2 in the master equation @ and 
using the normalisation condition Pqa + Paa = 1, we 
obtain the following equation. 



aP 



AA 



PP 2 aa + iPaa -6 = 



(5) 
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2beD + bD z 
D 2 + 
2beD + eP> 
2e 2 + 2eD. 



be 2 , 



2eD - AbeD - 2bD 2 ~ 2be 2 , 



e 2 b- 



2eW + 2WD, 



(G) 



Now Eq. © is in terms of two unkowns b and Paa and 
hence we must fix b in order to find Paa- It is known 
that in the limit D — > oo , the front moves with the Fisher 
velocity Vq = 2\feD. Also from equation 0), the front 
velocity in terms of p\ is given as V ~ Dpi = DPaa, 
when D is very large compared to e and W. Equating this 
velocity with the Fisher velocity Vo, we get Paa = 
Now substituting this value of Paa in Eq. JSJ, we have 
an equation which is linear in b and which in the limit 
D — > oo gives b = j. Substituting this value of b in the 
Eq. (J5J we obtain the following cubic equation in Paa- 



(c 2 



D 2 + 2eD)P\ A + (2e 2 + 2D 2 



AeD)Pi A + (5e 2 
+ 6eD + 8eW + 8WD)P AA ~ 8e 2 - 8eD = (7) 



One can easily solve the above cubic equation and the 
density at site just behind the front particle p\ = Paa 
can be obtained. The results obtained have been shown 
in the Figs. J3J and 10} and are in excellent agreement 
with the simulation results. Here we have also shown the 
percentage relative error in p\ i. e., Pl pS Pl ! x 100, where 
p\ and pi correspond to simulation and analytic results. 
Once we know this p\ , the front velocity is obtained from 
Eq. (JIJ and we have shown it in the Fig. (JSJ and |JBJ. 
Here, we also observe very good agreement with the sim- 
ulation results. In Figs. J3J andJSJ we have compared 
the results with that of Kerstein two particle represen- 
tation. The interesting thing that we notice here is that 



the results obtained from the present work is closer to 
the simulation results as compared to that obtained us- 
ing two particle representation. 
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FIG. 3: Comparision between simulation and analytic results 
for pi for e = 0.05, W = 0.0 and for different values of D. The 
open circle corresponds to Kerstein two particle representa- 
tion while the closed circle is the result from the present work. 
We note that the results of TSSCM are essentially coincident 
with the simulation results. Inset: Percentage relative error 
in pi. 
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FIG. 4: Comparision between simulation and analytic results 
for pi for e = 0.25, W = 0.125 and for different values of 
D. The closed triangle corresponds to the simulation results 
while the closed circle is the results from present work. We 
note that the analytic results are essentially coincident with 
the simulation results. Inset: Percentage relative error in pi. 

To conclude, we have developed a two site self consis- 
tent method for the propagating fronts in the reaction 
diffusion system A <-*• 2A whose results are in excellent 
agreement with the simulation results for all parameter 
regimes. We observe that for W — 0, the results ob- 
tained are better than that using Kerstein's two-particle 
representation. We also notice that the present work ap- 
pears to have an advantage over Kerstein's two-particle 
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representation due to two key factors : firstly, TSSCM 
doesnot neglect the spatial density correlation which, in- 
deed, is neglected in two particle representation by using 
the product measure ansatz , secondly, TSSCM forms 
a closed set of states under transitions due to the mi- 
croscopic processes for all parameter regimes while the 
two-particle representation does not provide a closed set 
for W 7^ 0. The simplicity of our analytic method pro- 
vides an scope to study the velocity of propagating front 
in other reaction diffusion processes as well. 



FIG. 5: Comparision between simulation and analytic results 
for front velocity V for e = 0.05, W = 0.0 and for different val- 
ues of D. The open circle corresponds to Kerstein two particle 
representation while the closed circle is the result from present 
work. We note that the results of TSSCM are essentially coin- 
cident with the simulation results. Inset: Percentage relative 
error in V. 
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FIG. 6: Comparision between simulation and analytic results 
for front velocity V for e = 0.25, W = 0.125 and for different 
values of D . The closed triangle corresponds to the sim- 
ulation results while the closed circle is the results from the 
present work. We note that the analytic results are essentially 
coincident with the simulation results. Inset: Percentage rel- 
ative error in V. 



[1] W. van Saarloos, Phys. Rep. 386, 29 (2003). 

[2] A. R. Kerstein J. Stat. Phys. 45, 921 (1986). 

[3] A. R. Kerstein , J. Stat. Phys. 53,703 (1988). 

[4] D. Panja, G. Tripathy and W. van Saarloos, Phys. Rev. 

E, 67, 046206 (2003) 
[5] Niraj Kumar and G. Tripathy, cond-mat / 0505599. 
[6] to be published. 

[7] R. A. Fisher, Ann. Eugenics , 7, 355 (1937). 



[8] E. Brunet and B. Derrida, Phys. Rev. E 56, 2597 (1997); 

J. Stat. Phys. 103, 269 (2001). 
[9] D. A. Kessler, Z. Ner and L. M. Sander, Phys. Rev. E 

58,107 (1998). 

[10] M. Bramson, Mem. Am. Math. Soc. 44, No. 285 (1983). 
[11] U. Ebert, W. van. Saarloos, Physica D 146, 1 (2000). 



